# Full text of "Unstable $m=1$ modes of counter-rotating Keplerian discs"

## See other formats

(N rn % Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 25 February 2013 (MN WT^ style file v2.2) Unstable m = 1 modes of counter— rotating Keplerian discs o. ^ ■ Mamta Gulati^'^'^, Tarun Deep Saini^''^ and S. Sridhar^'^ H ; C^ , ^ Indian Institute of Science, Bangalore 560 012, India ^ Raman Research Institute, Sadashivanagar, Bangalore 560 080, India ^ mgulati&rri.res.in , * tarun@physics.iisc. ernet. in , ^ ssridhar@rri.res.in 25 February 2013 o ^' ABSTRACT O . $—( , We study the linear m = 1 counter-rotating instability in a two-component, jrt ' nearly Keplerian disc. Our goal is to understand these slow modes in discs orbiting massive black holes in galactic nuclei. They are of interest not only be- cause they are of large spatial scale — and can hence dominate observations — but also because they can be growing modes that are readily excited by ac- > en , Cn . cretion events. Self-gravity being nonlocal, the eigenvalue problem results in a pair of coupled integral equations, which we derive for a two-component O ; softened gravity disc. We solve this integral eigenvalue problem numerically for various values of mass fraction in the counter- rotating component. The eigenvalues are in general complex, being real only in the absence of the counter-rotating component, or imaginary when both components have iden- tical surface density profiles. Our main results are as follows: (i) the pattern speed appears to be non negative, with the growth (or damping) rate be- ing larger for larger values of the pattern speed; (ii) for a given value of the pattern speed, the growth (or damping) rate increases as the mass in the counter-rotating component increases; (iii) the number of nodes of the eigenfunctions decreases with increasing pattern speed and growth rate. Ob- servations of lopsided brightness distributions would then be dominated by modes with the least number of nodes, which also possess the largest pattern speeds and growth rates. Key words: instabilities — stellar dynamics — celestial mechanics — galax- ies: nuclei © 0000 RAS 2 Gulati, Saini & Sridhar 1 INTRODUCTION Most galactic nuclei are thought to host massive black holes and dense clusters of stars whose structure and kinernatics are correlated with glo bal galaxy properties (JGebhardt et al. Ferrarese fc Merritt 2000 Gebhardt et al 1996 : 2000l ). Such correlations raise questions of great interest related to galaxy formation and the growth of nuclear black holes. The nearby large spiral galaxy M31 has a n off-centered peak in a double-peaked brightness distribution aroun d its nuclear black hole ( jLight et al. 1974 : Lauer et al.l 1993 1998 Kormendy fc Bender 19991 ). This lopsided brightness distribution could arise naturally if the a poapses o stellar orbits, orbiting the black hole, happened to be clustered together (JTremaine many 19951 ). Since then, kinematic and dynamical rnodels of such an eccentric disc have been con structed by several authors 2002 Peiris fc Tremaine Bacon et al. 2001 Salow &: Statler 2001 Sambhus fc Sridhai Sambhus fc Sridhar 20031 ). Of particular interest to this work is the model ofy (J2002l ). which included a few per cent of stars on retrograde (i.e. counter-rotating) orbits. They proposed that these stars could have been accreted to the centre of M31 in the form of a globular clu ster tha t spira led in due to dynamical friction. This proposal was motivated by the work of Toumal (120021 ) . who demonstrated that a Keplerian axisymmetric disc is susceptible to a linear lopsided instability in the m = 1 mode, even when a small fraction of the disc mass is in retrograde motion. Toumal (|2002| ) considered the linearized secular dynamics of particles orbiting a point mass, wherein particle orbits may be thought of as slowly deforming elliptical rings of small eccentricities. The m = 1 counter-rotating instability was studied analytically for a two-ring system, and numerically fo r a many-ring system. T he corresponding problem for continuous discs was then studied by Sridhar fc Saini (I2OIOI ). who proposed a simple model with dy- namics that could be studied largely analytically in the Wentzel-Kramers-Brillouin (WKB) approximation. Their model consisted of a two-component softe ned gravity d isc, orbiting a massive central black hole. Softened gravity was introduced by iMillen (jl97ll ) to simplify the analysis of the dynamics of stellar systems. In this form of interaction, the Newtonian 1/d gravitational potential is replaced by l/y/d'^ + 6^, where 6 > is called the softening length. In the context of waves in discs, it is well kn own that the softening len gth mimics the epicyclic radius of stars on nearly circular orbits ( iBinney &: Tremainell2008l ). Therefore, a disc composed of cold coUisionless matter, interacting via softened gravity, provides a surrogate for a hot coUisionless disc. © 0000 RAS, MNRAS 000, 000-000 Unstable m = 1 modes of counter-rotating Keplerian discs 3 Sridhar fc Saini (l2010f ) used a short-wavelength (WKB) approximation, derived analyt- ical expressions for the dispersion relation and showed that the frequency u is smaller than the Keplerian orbital frequency by a factor proportional to the small quantity e = Mct/M (which is the ratio of the disc mass to mass of the central object); in other words, the modes are slow. The WKB dispersion relation was used to argue that equal mass counter-rotating discs with the same surface density profiles (i.e. when there is not net rotation) could have unstable modes. They also argued tha t, for an arbitrary mass ratio, the discs must be unre- alistically hot to avoid an instability. ISridhar fc Saini I ( 120101 ) then used Bohr-Sommerfeld quantization to construct global modes, within the WKB approximation. A matter of con- cern is that the wavelengths of the modes could be of order the scale length of the discs; the modes being large-sca le it is possible t h at th e WKB approximation itself is invalid. Another limitation is that Sridhar fc Saini ( I2OIOI ) could construct (WKB) global modes only for the case of equal mass discs. Therefore it is necessary to address the full eigen- value problem to understand the systematic behaviour of eigenvalues and eigenf unctions. To this end, we formulate the eigenvalue problem for the linear, slow, m = 1 modes in a two-component, softened gravity, counter-rotating disc. Due to the long-range nature of gravitational interactions, we have to deal with a pair of coupled integral equations defining the eigenvalue problem. We draw some general conclusions and then proceed to solve the equations numerically for eigenvalues and eigenfunctions. In § 2 we introduce the unperturbed two-component nearly Keplerian disc, define the apse precession rates, discuss the potential theory for softened gravity, and derive the cou- pled, linear integral equations that determine the eigenvalue problems for slow m = 1 modes. A derivation of the relationship between the softened Laplace coefficients (used in the po- tential theory of § 2 and § 3) and the usual (unsoftened) Laplace coefficients is given in the Appendix. We specialize to discs with similar surface density profiles in § 3, when the two coupled equations can be cast as a single integral equation in a new mode variable; this allows us to draw some general conclusions about the eigenvalue problem. We also discuss in detail the numerical method to be employed. Our results are presented in § 4, where the properties of the stable, unstable and overstable modes are discussed. Conclusions are offered in § 5, where we seek to provide a global perspective on the correlations that occur between the pattern speeds, growth rates and eigenfunctions, as well as the variations of these quantities on the mass fraction in retrograde orbits. © 0000 RAS, MNRAS 000, 000-000 4 Gulati, Saini & Sridhar 2 FORMULATION OF THE LINEAR EIGENVALUE PROBLEM We consider linear non-axisymmetric perturbations in two counter-rotating discs orbiting a central point mass M. The discs are assumed to be coplanar and consist of cold collisionless particles which attract each other through softened gravity. However, the central mass and the disc particles interact via the usual (unsoftened) Newtonian gravity. Softened gravity is known to mimic the effects of velocity dispersion, so our discs are surrogates for hot stellar discs. We assume that the total mass in the discs, M^, is small in comparison to the central mass. Since e = Md/M ^ 1, the dynamics is dominated by the Keplerian attraction of the central mass, and the self-gravity of the discs is a small perturbation which enables slow modes. Below we formulate the linear eigenvalue problem of slow modes. 2.1 Unperturbed discs We use polar-coordinates r = (-R , 0) in the plane of the discs, with the origin at the loca- tion of the central mass. Throughout this paper the superscripts '-I-' and ' — ' refer to the prograde and retrograde components, respectively. The unperturbed discs are assumed to be axisymmetric with surface densities S^ (R). The disc particles orbit in circles with veloci- ties, v^ = ± RQ{R)e^, where f2(i?) > is the angular speed determined by the unperturbed gravitational potential. ^R) = -™ + UR) ■ (1) The first term on the right side is the Keplerian potential due to the central mass, and $d(-R) is the softened gravitational potential due to the combined self-gravities of both the discs: $,(r) = -G I Sl^2±S£2dV, (2) Ir — r'l + b"^ where b is the Miller softening length; the potential $d(-R) is 0(e) compared to GM/R. Test particles for nearly circular prograde orbits have azimuthal and radial frequencies, Q and k, given by © 0000 RAS, MNRAS 000, 000-000 Unstable m = 1 modes of counter-rotating Keplerian discs 5 n\R) K^R) GM 1 d^rf + R^ GM R di?' 3d$rf d=^$. (3) (4) i?3 ' i? di? di?2 ■ The line of apsides of a nearly circular eccentric particle orbit of angular frequency ±Q{R), subjected only to gravity, precesses at a rate given by ±zo{R), where w{R) = n{R) - k{R) 2 d RdR <l>d(i?) + Oie' (5) 2n{R) \RdR ' di?2 The cancellation of the 0(1) term, {GM/ R^), which is common to both Vl? and h? makes w ~ 0{e). This is the special feature of nearly Keplerian discs which is responsible for the existence of (slow) modes whose eigenfrequencies are ~ 0{e) when compared with orbital frequencies. 2.2 Perturbed discs Let vj(r,t) = ■u^(r,t)e/? + t>^(r,t)e0 and S^(r,t) be infinitesimal perturbations to the velocity fields and surface densities of the ± discs, respectively. These satisfy the following linearized Euler and continuity equations: -^ + (v,^ ■ V)v,± + (v± ■ V)v| -v$„, ^ + V ■ (S,^v,^ + S^,^ 0, (6) (7) where $a(i',^) is the perturbing potential. Fourier analyzing the perturbations in t and 0, we seek solutions of the form, Xa(r, t) = J2m -^Ti^) exp[i(m0 — ut)] . Then u mi: D± d {±mil — w)— — ± dR ,m± D± K^ d m i(±mf]-a;)S;^ © 0000 RAS, MNRAS 000, 000-000 m± 1 d RdR it„,»n±\ iR^,< 2mVL R <T)™ (8) -u) ^m (9) ir + 1 ^2± m± _ T d a = 0, (10) 6 Gulati, Saini & Sridhar where D, ± K^ — {±.mVt — oj) (11) The above equations determine m" r7i± „,m± and S^^ in terms of the perturbing potential $^ ; this would be the solution were $^ due to an external source. We are interested in modes for which $JJ^ arises from self gravity. In this case it depends on the total perturbed surface density, [SJ^^(i?) + S™~(_R)] . Manipulating the Poisson integral given in Eq. [21 we obtain K{R) R'dR'PUR,R') [sr(^') + sr (^')] (12) where the kernel PUR,R') =-^B[-\a,P) + "^^ (5m,l + Sm,-l) ■ (13) The second term on the right side is the indirect term arising from the fact that our coordi- nate system can be a non-inertial frame, because its origin is located on the central mass. The first term is the direct term coming from the perturbed self-gravity. Here i?< = min(i?, R') , i?> = max(_R, R') , a = R</ Ry and /3 = 6/-R> . The functions. i?f)(a,/3) =1 Tde- ^ Jo (1 cos m6 2apsl+ a2j^2)./2 ' (^^) are "softened Laplace coefficients", introduced in iToumal (l2002l ). They can be expressed in terms of the usual (unsoftened) Laplace coefficients, as shown in Appendix A. We note that the unperturbed disc potential ^^ can be obtained from the unperturbed disc density, S+(i?) + S^ (i?) , by using Eq. ([12]) with m = . 2.3 Slow m = 1 modes Modes with azimuthal wavenumber m = ±1 are slow in the sense that their eigenfrequencies, u, are smaller than the orbital frequency, Q, by a factor ^ 0(e) . Without loss of generality we may choose m = 1. In the slow mode approximation (JTremaind [200ll ) . we use the fact that fi > tj in Eqs. ([8])— ([U]), and write © 0000 RAS, MNRAS 000, 000-000 Unstable m = 1 modes of counter-rotating Keplerian discs 7 u i± ^ Df Y'dR ' d 2fi ,i± ± Df VL d Vt 2"Zr ^ 'r $1 $i where ±-^^K^ + ^^{RT^Wa^) + ^^^A^ Df = ±2n{uj T ro). 0, (15) (16) (17) (18) Eqs. flT^ and flT^ imply the following relations between the perturbed velocity amplitudes: u i± 2i^^ , m We use Eqs. ( !T9|) in the continuity equation (IT7|) to eliminate m^^ and write, ± OEi^ 2 d Combining Eqs. (112]), ([HD and (118])— (EO]) we obtain (R'^'^K^) ■ (20) [C^ T MR)]v'a^iR) °° di?'i?'i/2 r 5 2R^n{R') \dR d [R^Pi{R,R' dR , [i?'^/^Sl(i?>^(i?0 - i2'V2s^-(i?>i-(i?')] . (21) We rewrite this by defining zHR) n{R) 1/2 ,1± (22) use the fact that D.{R) oc R ^^^ for a Keplerian flow, and integrate by parts to obtain, © 0000 RAS, MNRAS 000, 000-000 8 Gulati, Saini & Sridhar /■OO 1 D/ [uj T vj{R)]z^{R) = - j ^2TiR,R') + / —2J^iR,R') j:j{R')j:^{Ry n{R')n{R) n{R')n{R) 1/2 z+{R 1/2 z-iR'), (23) where TiR, R!) 1 1 d l + l^^\Pi{R,B!) (24) 2d\nR'J \' ' 2d\nR^ It is convenient to write S^(i?) = ri{R)T,^{R) and S^(i?) = (1 — ri{R))T,^{R) , where r]{R) is the local mass fraction in the unperturbed counter- rotating component; by definition, ^ ?7(-R) ^ 1- Then, Eq. fl23|) can be recast as uz+{R) = +UJZ^{R) + dR R r[(i - ^(^'))(i - viR))]'^'}C{R, R')z^{R') dR R -[V{R'){1-V{R))Y/')C{R,R')Z~{R') UZ (i?) = -UJZ (R) + dR R -[{l-v{R')MR)]'^'}C{R,R')z^R') dR R yUR'UR)Y"lC{R,R')z-{R!) where the kernel (25) JC{R,R') 1/2 - 2 2nG ^AR')^<iiR) _ n{R')n{R) - ^AR')^ARy "" n{R')n{R) ^AR')^AR) n{R')n{R) J^{R, R') i.i- ' 2d\nR' l + l-J-^\PiiR,R') 2(9 In i? .1/2 ?(1). ' + 2 9h^J V + 2dh[R)^r^r~ • ^''^ Therefore the kernel /C(i?, i?') is a real symmetric function of R and i?'|j Using Eqs. ( IT9|) and ( l22l) . we can relate the eigenfunctions ;z'''(-R) and z~{R) to each other: ^ The contribution from the indirect term in Pi{R, R') vanishes, because (2 + d/d In R') R' ^ = . © 0000 RAS, MNRAS 000, 000-000 Unstable m = 1 modes of counter-rotating Keplerian discs 9 ^{l-r]{R))[u-w{R)]z+{R) = ^/^[u + w{R)]z-{R) . (27) This relation can, in principle, be used to eliminate one of z'^{R) or z^{R) from the coupled Eqs. (|25l) . in which case the eigenvalue problem can be formulated in terms of a single func- tion (which can be either z'^{R) or z~{R)). However, such a procedure results in a further complication: the eigenvalue, u, will then occur inside the R' integral in the combination, {u ±w)/{uj ^ w), and this makes further analysis difficult. Eqs. (!25|) are symmetric under the (simultaneous) transformations, {'+' ,ri{R) ,00} — )■ {'— ' , [1 — //(i?)] ,—uj}, which in- terchange the meanings of the terms prograde and retrograde. It seems difficult to obtain general results when E^(i?) and S^(i?) have different functional forms. Below we consider the case when the mass fraction, rj, is a constant; i.e. when both Sj'(i?) and S^(i?) have the same radial profile. 3 THE EIGENVALUE PROBLEM FOR CONSTANT r] DISCS When the counter-rotating discs have the same unperturbed surface density profiles, i.e. S^(i?) oc S^(i?) , some general results can be obtained. This case corresponds to the choice rj = constant, so that S^(_R) = rj'E^lR) and Sj"(i?) = {l — rj)'E^{R) . Then the eigenfunctions z~^{R) and z~{R) are related to each other by, [u-w{R)]^z+{R) = [uj + w{R)]^/{l-ri)z-{R). (28) Let us define a new function, Z{R), which is a linear combination of z^{R) and z~{R) : Z{R) = ^l-r]z+{R) - ^z-{R). (29) Then equations ( 12^ can be manipulated to derive a closed equation for Z{R) : w^ — w"^ /CO A T3I — ICiR,R')ZiR'), (30) [l-2r])u + ta_ We note that, in this integral eigenvalue problem for the single unknown function Z{R) , the (as yet undetermined) eigenvalue u occurs outside the integral. Once the problem is solved and Z{R) has been determined, we can use Eq. f l28|) and f l29|) to recover z^{R) : © 0000 RAS, MNRAS 000, 000-000 10 Gulati, Saini & Sridhar (1 — 277)0; + tJ7 (1 — 2?7)CJ + tJ7 Some general conclusions can be drawn: (i) In Eq. ( |30ll . the kernel )C{R,R') is real symmetric. Therefore, the eigenvalues, u, are either real or come in complex conjugate pairs. (ii) When n = 0, the counter-rotating component is absent, which is the case studied by Tremaind (1200 if ): then the left side of Eq. (!30|) becomes {uj — zo) Z . Since the kernel 1C{R, R') is real symmetric, the eigenvalues uj are real, so the slow modes are stable and oscillatory in time. Then the eigenfunctions, Z{R) may be taken to be real. Therefore z^{R) = Z{R) is a real function, and z~(R) = . g (iii) When rj = 1/2, there is equal mass in the counter-rotating component, and the surface densities of the ± discs are identical to each other. This case may also be thought of as one in which there is no net rotation at any radius. Then Eq. (|30|) becomes w-\R){u^-w\R))Z{R) = / ^}C{R,R')Z{R'), (32) Jo R Since the kernel /C(i?, R') is real symmetric, w^ must be real. There are two cases to consider, when the eigenvalues, w, are either real or purely imaginary. • When UJ is real, the slow modes are stable and oscillatory in time. The eigenfunctions z'^{R) can be taken to be real functions. • When u is imaginary, the eigenvalues come in pairs that are complex conjugates of each other, corresponding to non-oscillatory growing/damped modes. Let us set r] = 1/2 and u = i'j (where 7 is real) in Eq. (ISTl) : The function Z{R), which is a solution of Eq. ( l32|) . can be taken to be a real function multiplied by an arbitrary complex constant. It is useful to note two special cases: (i) when Z{R) is purely imaginary, then z~^{R) and z~{R) are complex conjugates of each other; (ii) when Z{R) is real, z~^{R) is equal to minus one times the complex conjugate of z~{R). To make progress for other values of r^, it seems necessary to address the eigenvalue problem numerically; the rest of this paper is devoted to this. ■^ When 'rj = 1, the eigenvalues, oj, are again real, with z^ (R) = Z{R) a real function, and z'^{R) = . © 0000 RAS, MNRAS 000, 000-000 Unstable m = 1 modes of counter-rotating Keplerian discs 11 3.1 Application to Kuzmin discs For numerical explorations of the eigenvalue problem, we consider the case when both the unperturbed ± discs are Kuzmin discs, with similar surface density profiles: S^(_R) = rjT,^{R) and S+(i?) = (1 - r])T.^{R) , where ^AR) aM, (34) 27r(i?2 + a2)3/2 is the total surface density, Md is the total disc mass and a is the disc scale length. Kuzmin discs, being centrally concentrated, are re asonable candidates for unperturbed dis cs. More- over, earlier investigations of slow modes ( iTremaine 2001 Sridhar fc Saini 2010h have ex- plored modes in Kuzmin discs, so we find this choice useful for comparisons with ear- lier work. The characteristic values of orbital frequency and surface density are given by VL* = ^GM ja^^ and S^ = Md/a^, respectively. The coupled Eqs. ( I25l) can be cast in a di- mensionless form in terms of these physical scales. The net effect is to rescale the eigenvalue u to a, where a Q*a u (35) In the following section, all quantities are to be taken as dimensionless; however, with some abuse of notation, we shall continue to use the same symbols for them. 3.2 Numerical method Our method is broadly similar to iTremaind (1200 ll ). In order to calculate eigenvalues and eigenf unctions numerically, we approximate the integrals by a discrete sum using an Ap- point quadrature rule. The presence of the term dR/i? in the integrals suggests that a natural choice of variables is m = log(-R) and v = log(-R'), where, as mentioned above, R stands for the dimensionless length R/a. The use of a logarithmic scale is numerically more efficient, because it induces spacing in the coordinate space that increases with the radius. This handles naturally a certain expected behaviour of the eigenf unctions: since the surface density in a Kuzmin disc is a rapidly decreasing function of the radial distance, we expect the eigenfunctions to also decrease rapidly with increasing radius. Therefore, discretization of the coupled equations fl25|) follows the schema: © 0000 RAS, MNRAS 000, 000-000 12 Gulati, Saini & Sridhar (36) where Wj's are suitably chosen weights. Then the discretized equations can be written as a matrix eigenvalue problem: where AC = ^C: (37) 7])wjlCij + Wj5i_ ■\/ri{l -vi)wjlCi ^/riiY^^WjlCi -rjWjK,, «j Wjbij and V z- z, (38) The 2N X 2N matrix A has been represented above in a 2 x 2 block form, where each of the 4 blocks is a A^ x A^ matrix, with row and column indices i and j. Note that 5ij is the Kronecker delta symbol, and no summation is implied over the repeated j indices. Thus we have an eigenvalue problem for eigenvalues a, and eigenvectors given by the 2N dimensional column vector C • The use of unequal weights destroys the natural symmetry o f the kernel, but t lis is readily restored through a simple transformation given in § 8.1 of Press et al. (Il992[ ). The grid for our numerical calculations covers the range —7 ^ logi? ^ 5, which is divided into A^ = 4000 points; large r values of N give similar results. We note some differences with Tremaind (l200l[ ) concerning details of the numerical Tremaine metho d and assumptions. The major difference is in the treatment of softening: In teOOll ). a dimensionless softening parameter /3 = h/R was introduced, and the eigenvalue problem for slow modes was solved by holding the parameter /3 constant. This renders the physical softening length, 6, effectively dependent on radius, making it larger at larger radii, thereby not corresponding to any simple force law between two disc particles located at different radii. We have preferred to keep h constant, so that the force law between two disc parti cles is through th e usual Miller prescription. Other minor differences in treatment are: 1 m Tremaind ( 120011 ). the disc interior to an inner cut-off radius was assumed to be frozen. In contras t we u se a straightforward inner cut-off radius of 10 ^ , as mentioned above; (ii) Tremaind ( 120011 ) uses a uniform grid in logi?, with four-point quadrature in the intervals © 0000 RAS, MNRAS 000, 000-000 Unstable m = 1 modes of counter-rotating Keplerian discs 13 c O o 1 — I I I lllll 1 — I I I Mill 1 — I I I Mill 1 — I I I illlj 1 — I I I I I — I I I illlj 1 — I I illlll 1 — I I I Mill 1 — I I I llllj 1 — TTT IT] 1 — I I I llll| 1 — I I I llil| 1 — I I I iill| a = -0.2370 I I m ill — I I m ill — I I mm — i i mm — h-h = -0.1820 ] 1 I I I llll| 1 — I I I llll| 1 — I I I llll| CT = -0.1422 I I lllllll — I I lllllll — I I I Illlll — I I lllllll — H-H __ CT = -0.1101 » 0.01 0.1 1 10 0.01 0.1 1 10 Radial coordinate R Figure 1. Slow g-modes in a single, prograde (77 = 0) Kuzniin disc with A = 0.1 and b = 10~ . The panels are labeled by the scaled eigenvalue cr. between consecutive grid points; we also use a uniform grid in logi? but instead employ a single A^-point quadrature for integration. Were we dealing with unsoftened gravity (i.e. the case when 6 = 0), the diagonal elements of the kernel would be singular. Hence, when the softening parameter b is much smaller than the grid size, accuracy is seriously compromised by round-off errors. Typically, the usable lower limit for b is ~ 10 ~^. 4 NUMERICAL RESULTS We obtained the e igenvalues and eigen: 'unctions of equation ( 1371) using the linear algebra package LAPACK (jAnderson et al.lll999h . We now present the results of our calculations for specific values of rj. As noted earlier, interchanging the meaning of prograde and retrograde © 0000 RAS, MNRAS 000, 000-000 14 Gulati, Saini & Sridhar b = 0.05 h a = 1.0357 0) -1 -2 c a; o u ■1 - -2 I I I lllll| 1 I I lllll| 1 I I I llll|| 1 I I lllll| 1 — I I I I I II lllll| 1 I I lllll| 1 I I I llll| 1 I I lllll| 1 — TT I I m ill — I I m ill — I I mm ' i i mm — h44+ b = 0.1 + a = 0.4181 b = 0.2 h a = 0.1775 uL uL -AJ I I lllllll — I I lllllll — I I lllllll ' I I lllllll — H4+ b = 0.4 + a = 0.1147 uL iJ_ 0.01 0.1 10 0.01 0.1 10 Radial coordinate R Figure 2. Slow p-modes in single, prograde (77 = 0) Kuzniin disc with no external source (i.e A = 1). The panels are labeled by the scaled eigenvalue a, and the softening parameter, b (which has been scaled with respect to a, the disc scale length). orbits leave the results invariant under the transformation (^7,0;) — )■ (1 — r/, — w); therefore, we present results below only for ^ r/ ^ 1/2 . 4.1 No counter-rotation: 77 = es rotate i n the prograde sense. The eigen- Tremaind (1200 ll ). who also showed that the We are dealing with a single disc whose partic value problem for this case was studied first by eigenvalues are real; in other words, the disc supports stable slow modes. We consider this case first to benchmark our numerical method as well as assess the differences in results that may arise due to the mann er in which softe ning is treated. To facilitate comparison we use the same nomenclature as Tremaind (J200lh . Briefly, modes corresponding to positive and negative eigenvalues are referred to as "p-modes" and "g-modes" , respectively; we also © 0000 RAS, MNRAS 000, 000-000 O : CD -M D q - d : O to I O Unstable m = 1 modes of counter-rotating Keplerian discs 15 T — I — I I I I I 1 — I I I 1 1 -i — I I I 1 1. = 0.05: J I J 10 100 1000 Number of nodes n Figure 3. Growth rate versus number of nodes for r) = 0.5, for two values of softening, b = 0.1 and b = 0.05. introduce a parameter A = (1 + /)^^, where / is a constant that mimics additional precession due to an external source of the form tJ7e(_R) = fwdi^R)] we define eccentricity e^ through ex GM -1/2 (39) and use the normalization, dR eliR) 1 (40) Our results corresponding to g-modes, for A = 0.1 and /3 = 10~^, are presented in Fig. [H where we plot modes with three or fewer nodes and give their eigenvalues. Results for p-modes for A = 1 (no external source) and various values of softening param eter b, are displayed in Fig. [2l These figures are to be compared with Fig. 3 and Fig. 6 of © 0000 RAS, MNRAS 000, 000-000 Tremaine 16 Gulati, Saini & Sridhar 2 1 Lb = 0.05 1 ^1 -1 -2 1 -1 1 — I — I — I JL b = 0.1 I I I iiiii| — I I I ii mL — I I I iiiii| — I I I iiiii| b = 0.05 V T = 0.0282 I I I iiiii| — I I I iiiii ]| I I I iiiii| — I I I iiiii| d I rri I I I = 0.0355 I I I iiiii| — I I 1 1 1 1 I 'l l — I I I iiiii| — I I I iiiii| uJ I 0.1 1 10 100 0.1 1 Radial coordinate R 10 100 Figure 4. Eigcnfunctions Z{R) are plotted as a function of the radial coordinate R, for rj = 0.5 . The panels are labeled by the values of the growth rate, T, and softening, h. Figure 5. Gray-scale plots of surface density perturbations, E^ (R,<j>,t) at time t = for the parameter values, r] = 0.5 and fe = 0.1 , and r = 0.0230. White/black correspond to the maximum positive/negative values of the perturbations. © 0000 RAS, MNRAS 000, 000-000 Unstable m = 1 modes of counter-rotating Keplerian discs 17 m (120011): the e j genf u nctions are of broadly similar form, but the eigenvalues differ from those Tremaind (120011 ) by upto ~ 30% . 4.2 Equal counter-rotation (or no net rotation): rj = 1/2 This case was studied by Sridhar fc Saini (]2010f ). who derived the following local or "WKB" dispersion relation: U) zu w n{R) k\ exp{—\k\b) (41) From this expression they concluded: if zu happens to be positive then u is real (and the disc is stable), but zu is negativ e for most continuous d iscs which implies that u can be either real or purely imaginary. Sridhar fc Saini ( I2OIOI ) also studied global modes using Bohr-Sommerfeld quantization, which will be discussed later in this section. We have proved in the last section that the eigenvalues are either real (stable oscillatory modes), or purely imaginary (non-oscillatory growing and damped modes). Here we focus on the growing (unstable) modes; we define the growth rate of perturbations as. in order to facilitate comparisons with b\(T\ (42) Sridhar fc Saini mm. In Fig. IS l we p lot F for b Sridhar &: Saini torn found two equal to 0.1 and 0.05, versus the number of nodes of Z(R). separate branches in the spectrum, corresponding to long and short wavelengths (for each value of b). Comparing our Fig. |3]with their Figs. 3 & 4, we see that our results are more consistent with their short-wavelength branch than with their long-wavelength branch. This disagreement is probably because the long-wavelength branch corresponds to kR ~ 1, where WKB approximation breaks down. Moreover, the agreement between our results and their short-wavelength branch holds only in a broad sense, because there are differences in the numer ical values of the eigenvalues. We trace this difference to the fact that ISridhar fc Saini ( 120101 ) used an analytical result for the precession frequency corresponding to unsoftened gravity, whereas we have consistently used softened gravity for all gravitational interactions between disc particles. This probably also results in another difference between our results: according to Sridhar fc Saini ( I2OIOI ). 10"'^ < F < 10~^; however, as can be seen from our Fig. [31 we obtain values of F both inside and outside this range. Changing the value of b causes a horizontal shift in the spectrum, which is consistent with their results. © 0000 RAS, MNRAS 000, 000-000 18 Gulati, Saini & Sridhar a 0.2 0.1 - -0.1 -0.2 0.1 -0.1 -0.2 0.1 -0.1 h -0.2 "1 ' r -n = 0.1 -^r 0.25 4- H I . I - I T] = 0.4 _L -0.2 0.2 0.02 -0.02 0.02 b" -0.02 0.02 -0.02 —1 1 1 1 1 1 1 1 r- 7] = 0.1 -^S: H 1 1 h 0.25 H 1 1 1- 4 H 1 1 h 0.4 -+ — I — I — h _I I I L_ -0.05 0.05 Figure 6. Distribution of eigenvalues in the complex a plane, for r] = 0.1,0.25 and 0.4. Panels labeled a give an overview, whereas the panels labeled b provide an close-up view near the origin. In Fig m we plot a few of these eigenfunctions as a function of R for both values of b equal to 0.1 and 0.05. Note that, from the discussion in the previous section, Z{R) can always be chosen to be a real function of R. The smallest number of nodes corresponds to the largest value of the growth factor. The eigenfunctions with the fewest nodes have significant amplitudes in a small range of radii around -R ~ 1, and this range increases with the number of nodes (and correspondingly, the growth rate decreases). Figgis a gray-scale plot of the surface density perturbations in the ± discs, S^(/2, 0, t = 0). Note the relative phase shift between the it perturbations. For other value of b shown in Fig H] we get similar patterns. © 0000 RAS, MNRAS 000, 000-000 T 1 I I I I III 1 1 I I I INI 1 1 I I I INI 1 1 I I I I 1 11 2 F 1 N ■1 - -2 t n — I II Mill 1 - N3 -1 -2 Unstable m = 1 modes of counter-rotating Keplerian discs 7) = 0.1 19 £7^=0.2501 CTj=0.1754 n — I I I I Mil 1 — I I I I Mil 1 — I I I 1 1 III 1 — I I I 1 1 1 11 CTr=0.2221 CTj=0.1524 ^ I I I MM 1 1 I I I I 1 11 CTr=0.1574 CTj=0.2635 -q = 0.25 n — I I 1 1 MM I I I 1 1 III 1 — I I I I II 11 (Tr=0.1408 aj=0.2332 n 1 I I I I III 1 1 I I I MM 1 1 I I I MM 1 1 I I I I 1 11 77 = 0.4 CTp=0.0631 CTj=0.3011 r n — I I 1 1 MM — I I 1 1 1 ii| 1 — I I 1 1 II ij. (7(^=0.0565 (Tj=0.2670 ^ qL 0.1 10 100 0.1 Radial coordinate R 10 100 Figure 7. Real parts of the "most unstable" eigenfunctions Z{R), plotted as a function of the radial coordinate, R for b = 0.1 and for r] = 0.1, 0.25 and 0.4. Panels are labeled by the real and imaginary parts of the eigenvalues. 4.3 Other values of r] We present results for values of rj other t han and 1/2. T only because they were not explored by Sridhar fc Saini l ese ar e particularly interesting, not (J2OIOD, but because the eigenvalues can be truly complex, corresponding to growing and damped modes which precess with steady pattern speeds. We write the eigenvalues as a = aR + io"/. In Fig|6], we display the eigenvalues in the complex a-plane, for softening parameter 6 = 0.1 and for rj equal to 0.1, 0.25 and 0.4. Panels on the left, labeled (a), provide an overview, whereas the panels on the ues n ear the origin right, labeled (6), provide a close-up view of the distribution o f eigenva of the complex a-plane; this distribution is similar to Fig. 3 of iToumal ( 20021 ). We are able © 0000 RAS, MNRAS 000, 000-000 20 Gulati. Saini & Sridhar Figure 8. Gray-scale plots of surface density perturbations, E^ (_R, <f>, t) at time t = for the parameter values, r] = 0.25 and b = 0.1 , and cr = 0.1408 + iO.2332. White/black correspond to the maximum positive/negative values of the perturbations. ^^ o CTr = 1.263x10"^ CTj = 1.267x10"^ - A A 11 , , 1 1 i i i i 1 C7„ = 1.256x10"^ (7i = 3.006x10"^ 1 1 10 Radial coordinate R (7r = 5.723x10"^ (7r = 5.727x10"^ . (T, = 5.35x10"^ 11 1 0-, = 7.58x10"* . 11 r III 10 Radial coordinate R Figure 9. Real parts of two pairs of eigenfunctions Z{R) (from two arms of a branch), plotted as a function of the radial coordinate R, for b = 0.1 and i) = 0.25. Panels are labeled by the real and imaginary parts of the eigenvalues. © 0000 RAS, MNRAS 000, 000-000 Unstable m = 1 modes of counter-rotating Keplerian discs 21 to provide much more detail, essentially because we are dealing with continuous discs rather than a finite number of rings. As 77 incre ases from 0, t he eigenvalues go from real to complex, a bifurcation that has been traced in iToumal ( 12002| ) to a phenomenon identified by M. J. Krein due to the resonant crossing of stable modes. The complex eigenvalues come in complex conjugate pairs, so there are two branches to the distribution. As rj increases, the branches progressively separate and, for ?7 = 1/2 must lie along the positive and negative imaginary axes. It is intriguing that each of these two branches consists of more than one arm. In the close-up views provided by the (b) panels, it appears as if each of the branches has two arms; however, more detailed investigations are required to determine if there are more arms. The arms of each of the branches are most widely separated when rj = 0.25, which is the value of rj exactly midway in its range ^ 77 ^ 0.5 . The separations decrease as rj approaches either or 1/2; this is natural because, for rj = both branches must lie on the real axis and, for rj = 1/2 both branches must lie on the imaginary axis. The eigenfunctions are in general complex, and have a rich structure as functions of their eigenvalues. Since our interest is in the unstable modes, we now display in Fig. [7| plots of the Zji = ^[Z{R)] in Eq. fl29l) . corresponding to the "most unstable" modes (for softening parameter b = 0.1 and for t] equal to 0.1, 0.25 and 0.4). In other words, for some chosen value of aR, we display the real part of the eigenf unction corresponding to the largest value of aj. For a fixed value of 77, the number of nodes of the eigenfunctions decreases with increasing pattern speed and growth rate. Fig [8] is a gray-scale plot of the surface density perturbations in the ± discs, S^(i?, 0, t = 0), for the parameter values rj = 0.25 and 6 = 0.1, and a = 0.1408 + 10.2332 . It is also of interest to ask how eigenfunctions from two different arms of the same branch behave. To do this, we picked two eigenfunctions with nearly the same value of a^, but with values of aj corresponding to two different arms of one branch; Fig. [9] shows two such pairs of eigenfunctions for 77 = 0.25. We have looked at pairs of such eigenfunctions for other values of ?7, but do not display them, here we note what seems to be a general trend: (i) the two members of a pair are more similar to each other when the values of their an are closer to each other; (ii) the member of a pair with the smaller value of ai is more displaced toward larger radii. © 0000 RAS, MNRAS 000, 000-000 22 Gulati, Saini & Sridhar 5 CONCLUSIONS We study linear, slow m = 1 modes in softened gravity, counter-rotating Keplerian discs. The eigenvalue problem is formulated as a pair of coupled integral equations for the ± modes. We then specialize to the case when the two discs have similar surface density pro- files but different ± disc masses. It is of great interest to study the properties of the modes as a function of rj, which is the fraction of the total disc mass in the retrograde population. Recasting the coupled equations as a single equation in a new modal variable, we are able to demonstrate some general properties: for instance, when rj = 1/2, the eigenvalues must be purely imaginary or, equivalently, the modes are purely unstable. In other words, when the ± discs have identical surface density profiles then t here are grow i ng m = 1 modes with zero p a ttern speed, a conclu s ion w h ich is consistent with fll990h : Sellwood fc Merritt Araki](ll987f): fll994[ ): iLovelace et all (119971 1: iToumal ( 120021 ): iTremaind (12005f ). Palmer fc Papaloizou To study modes for general values of rj, the eigenva l ue pr oblem needs to be solved numer- ically. Our method is broadly based on iTremaind ( 120011 ). but there are some differences whose details have been discussed i n the text . The main point of departure is in the way that softening has been treated. In Tremaind (1200 ll ). a dimensionless softening parameter /3 = b/R was introduced, and the eigenvalue problem was solved by holding the parameter P constant. This procedure renders the physical softening length, b, effectively dependent on radius (making it larger at larger radii), thereby not corresponding to any simple force law between two disc particles located at different radii. We have preferred to keep b constant, so that the force law between two disc particles is through the usual Miller prescription. We calculate eigenvalues and eigenf unctions numerically for discs with surface density profiles of Kuzmin form. Kuzmin discs, being centrally concentrated, are r easonable candi dates for unperturbed discs. Moreover, earlier investigations of slow modes (JTremaine Sridhar fc Saini 2001 : 2010l ) have explored modes in Kuzmin discs, so this choice is p articularly usefu Tremaine for comparisons with earlier work. Comparing our results with those of ( 1200 ll ) for 1] = (when the slow modes are stable), we find that the eigenfunctions are of broadly similar form, but the eigenvalues differ by up to ~ 30% ; this is a result of the differ- ent ways in which we have treated softening. For the case of no net rotation (77 = 1/2), we find that the growth rates (of the unstable modes) we calculate are broadly co nsistent with the Sridhar fc Saini short- wavelength branch of the global WKB modes determined earlier by (I2OIOI ). but not their long-wavelength branch. This disagreement probably arises because © 0000 RAS, MNRAS 000, 000-000 Unstable m = 1 modes of counter-rotating Keplerian discs 23 the long- wavelength branch corresponds to wavelengths of order the disc scale length where WKB approximation breaks down. Moreover, the agreement between our results and their short-wavelength branch holds only in a broad sense, because there are differences in the numerical values of the eigenvalues. We trace this difference to the fact that Sridhar fc Saini ( I2OIOI ) used an analytical result for the precession frequency corresponding to unsoftened gravity, whereas we have consistently used softened gravity for all gravitational interactions between disc particles. We have also investigate eigenmodes for values of t] other than and 1/2. These cases ar e Sridhar &: Saini mm, particularly interesting, not only because they were not explored by but because the eigenvalues can be truly complex, corresponding to growing (and damped) modes with non zero pattern speeds. We have presented results for rj = 0.1,0.25 and 0.4 in the previous sections. Based on these, we interpolate and offer the following conclusions about the properties of the eigenmodes and their physical implications, for all values of r) (which is the mass fraction in the retrograde population): (i) For a general value of 77 (between and 1/2), the distribution of eigenvalues in the complex plane has two branches. These branches are symmetrically placed about the real axis, because the eigenvalues come in complex conjugate pairs. (ii) The pattern speed appears to be non negative for all values of rj, with the growth (or damping) rate being larger for larger values of the pattern speed. (iii) For a fixed value of rj, the number of nodes of the eigenfunctions decreases with increasing pattern speed and growth (or damping) rate. (iv) For a value of pattern speed in a chosen narrow interval, the growth (or damping) rate increases as rj increases from to 1/2. (v) Each of the two branches in the complex (eigenvalue) plane has at least two arms. When rj = 0, the eigenvalues are all real, so both branches lie on the real axis, with zero spacing between the arms. As 77 increases, the branches lift out of the real axis, and the arms separate. It appears as if the maximum separation between the arms happens when rj = 1/4. As rj increases further, the branches continue to rise with greater slope, while the arm separation begins decreasing. Finally, when 77 = 1/2, the arm separation decreases to zero as the branches lie on the imaginary axis. Observations of lopsided brightness distributions around massive black holes are some- what more likely to favour the detection of modes with fewer nodes than modes with a large © 0000 RAS, MNRAS 000, 000-000 24 Gulati, Saini & Sridhar number of nodes, because the former suffer less cancellation due to finite angular resolution. From items (ii) and (iii) above, we note that the modes with a small number of nodes also happen to be those with larger values of the pattern speed and growth rate, both qual- ities that enable detection to a greater degree. Having said this, it would be appropriate to note some limitations of our work. Softened gravity discs are, after all, surrogates for discs composed of collisionless particles (such as stars) with non zero thickness and velocity dispersions. It is necessary to formulate the eigenvalue problem for truly collisionless discs, in order to really deal with stellar discs around massive black holes. Meanwhile, our results will serve as a benchmark for future investigations of modes in these more realistic models. REFERENCES Anderson et. al., 1999, Society for Industrial and Applied Mathematics, 3rd ed., Araki, S. 1987, Astron. J., 94, 99 Bacon, R., Emsellem, E., Combes, F., Copin, Y., Monnet, C, & Martin, P. 2001, Astron. Astrophys., 371, 409 Binney, J., and Tremaine, S. 2008, Galactic Dynamics (2ed., Princeton: Princeton Univer- sity Press) Ferrarese, L., & Merritt, D. 2000, Astrophysical. J. Letters, 539, L9 Gebhardt, K., et al. 1996, Astron. J., 112, 105 Gebhardt, K., et al. 2000, Astrophysical. J. Letters, 539, L13 Kormendy, J., & Bender, R. 1999, Astroph. J. , 522, 772 Lauer, T. R., et al. 1993, Astron. J., 106, 1436 Lauer, T. R., Faber, S. M., Ajhar, E. A., Grillmair, C. J., & Scowen, P. A. 1998, Astron. J., 116, 2263 Light, E. S., Danielson, R. E., & Schwarzschild, M. 1974, Astroph. J. , 194, 257 Lovelace, R. V. E., Jore, K. P., & Haynes, M. P. 1997, Astroph. J. , 475, 83 Miller, R. H. 1971, Astroph. Space Science, 14, 73 Murray, C. D., and Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press) Palmer, P. L., & Papaloizou, J. 1990, Mon. Not. Roy. Ast. Soc, 243, 263 Peiris, H. V., & Tremaine, S. 2003, Astroph. J. , 599, 237 © 0000 RAS, MNRAS 000, 000-000 Unstable m = 1 modes of counter-rotating Keplerian discs 25 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Cambridge: University Press, — cl992, 2nd ed., Salow, R. M., & Statler, T. S. 2001, Astrophysical. J. Letters, 551, L49 Sambhus, N., & Sridhar, S. 2002, Astron. Astrophys., 388, 766 Sellwood, J. A., & Merritt, D. 1994, Astroph. J. , 425, 530 Sridhar, S., & Saini, T.D. 2010, Mon. Not. Roy. Ast. Soc, 404, 527 Touma, J. R. 2002, Mon. Not. Roy. Ast. Soc, 333, 583 Tremaine, S. 1995, Astron. J., 110, 628 Tremaine, S. 2001, Astron. J., 121, 1776 Tremaine, S. 2005, Astroph. J. , 625, 143 APPENDIX A: EXPRESSING SOFTENED LAPLACE COEFFICIENTS IN TERMS OF (UNSOFTENED) LAPLACE COEFFICIENTS Softened Laplace coefficients were defined in Touma (120021 ) as where We now write One solution for 7 and 6 is. 2 r B^{a,P) = - de TT Jo COS m9 A = l + a^ + P^ -2acos9 A = 7" + (5' -275 cos ^ (Al) (A2) (A3) 7 l + a^ + Z^^ , 1 /-— — „ , ^„,^ -—, + - V(l + a^ + /32)2 - 4a2 1/2 6 = ^ 7 Therefore (A4) i?r(«,/3) = l-'bT/2iSh) (A5) © 0000 RAS, MNRAS 000, 000-000 26 Gulati, Saini & Sridhar where &r/2(«) - I de- ll cos mO :i + «2 2a cos 0) s/2 ' a < 1 (A6) are the famihar (unsoftened) Laplace coefficients (jMurray &: Derniottlll999l ). From eqn. (1A5I1 . we must have ((5/7) < 1. That this is indeed true can be proved using eqns. (JA4]): 7 is a monotonically increasing function of /3^, hence 7^1, and (5/7) = (a/7^) < 1 . © 0000 RAS, MNRAS 000, 000-000